rm(list = ls())

##

library(tidyverse)
library(sf)

## Get treatment

treat <- read_rds('data/Data_Clean.rds') %>%
  dplyr::select(county_id, treat_categ) %>% 
  distinct(county_id, .keep_all = T) %>% 
  mutate(treat_categ = ifelse(is.na(treat_categ), 'none', treat_categ))

## Recode treatment

treat <- treat %>% 
 mutate(treat_categ = dplyr::recode(treat_categ, 
                                     `sehr schwer` = 'Highly affected',
                                     `schwer` = 'Highly affected',
                                     `schwach` = 'Weakly affected',
                                     `none` = 'Unaffected')) %>% 
  mutate(treat_categ = factor(treat_categ, levels = unique(treat_categ)[c(1,2,3)]))

#### Figure 2 - map #### 

## Get county shapefile

fname <- 'data/shapefiles/vg2500_krs.shp'
shp <- sf::st_read(fname)

## Merge

shp <- shp %>% left_join(treat, by =c('RS' = 'county_id')) %>% 
  filter(substr(RS, 1, 2) %in% c('05', '07'))

## Get state shapefile

states <- sf::st_read('data/shapefiles/vg2500_bld.shp')

## Reorder

shp <- shp %>% 
  mutate(treat_categ = factor(treat_categ, 
                              levels = unique(treat_categ)[c(2,1,3)]))

## Map

p <- ggplot(shp) +
  # municipality polygons
  geom_sf(data = shp, 
          aes(fill = treat_categ, color = treat_categ)) +
  geom_sf(data = states ,
          fill = NA, color = 'grey60') +
  geom_sf(data = shp %>% filter(RS == 	'07131'),
          fill = NA, color = 'white') +
  theme_minimal() +
  scale_fill_grey(name = 'Flooding intensity', ) +
  scale_color_grey(name = 'Flooding intensity') +
  theme(legend.position = 'bottom') +
  guides(fill=guide_legend(ncol=2))
p
